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[i]  Hydraulic  tomography  is  a  powerful  technique  for  characterizing  heterogeneous 
hydrogeologic  parameters.  An  explicit  trade-off  between  characterization  based  on 
measurement  misfit  and  subjective  characterization  using  prior  information  is  presented. 

We  apply  a  Bayesian  geostatistical  inverse  approach  that  is  well  suited  to  accommodate  a 
flexible  model  with  the  level  of  complexity  driven  by  the  data  and  explicitly  considering 
uncertainty.  Prior  information  is  incorporated  through  the  selection  of  a  parameter 
covariance  model  characterizing  continuity  and  providing  stability.  Often,  discontinuities 
in  the  parameter  field,  typically  caused  by  geologic  contacts  between  contrasting  lithologic 
units,  necessitate  subdivision  into  zones  across  which  there  is  no  correlation  among 
hydraulic  parameters.  We  propose  an  interactive  protocol  in  which  zonation  candidates  are 
implied  from  the  data  and  are  evaluated  using  cross  validation  and  expert  knowledge. 

Uncertainty  introduced  by  limited  knowledge  of  dynamic  regional  conditions  is  mitigated 
by  using  drawdown  rather  than  native  head  values.  An  adjoint  state  formulation  of 
MODFLOW-2000  is  used  to  calculate  sensitivities  which  are  used  both  for  the  solution  to 
the  inverse  problem  and  to  guide  protocol  decisions.  The  protocol  is  tested  using  synthetic 
two-dimensional  steady  state  examples  in  which  the  wells  are  located  at  the  edge  of  the 
region  of  interest. 

Citation:  Fienen,  M.  N.,  T.  Clemo,  and  P.  K.  Kitanidis  (2008),  An  interactive  Bayesian  geostatistical  inverse  protocol  for  hydraulic 
tomography,  Water  Resour.  Res.,  44,  W00B01,  doi:10.1029/2007WR006730. 


1.  Introduction  and  Background 

[2]  A  principal  challenge  in  hydrogeology  is  the  need  to 
obtain  details  on  characteristics  and  properties  of  an  aquifer 
over  potentially  large  areas.  The  installation  of  wells  is 
expensive  and  often  constrained  by  surface  constmction, 
property  access  considerations,  and  a  host  of  other  issues. 
Best  approaches  are  those  that  make  the  most  of  the  data 
available  within  the  context  of  errors.  Early  aquifer  charac¬ 
terization  was  based  on  pumping  tests  performed  in  a  single 
well  (converging  radial  flow),  assuming  a  homogeneous 
aquifer,  and  with  various  other  assumptions,  fitting  an 
analytical  solution  to  the  drawdown  curve.  This  work  began 
with  Slichter  [1899]  and  Thiem  [1906]  for  steady  state  and 
Theis  [1935]  extended  it  to  transient  conditions.  Compendia 
of  special  cases  under  various  conditions  are  presented  in 
textbooks  [e.g.,  Batu ,  1998;  Kruseman  and  De  Bidder, 
1990].  These  solutions  provide  a  homogeneous  bulk  aver¬ 
age  of  the  target  hydraulic  parameter.  However,  no  medium 
is  tmly  homogeneous  and  even  subtle  heterogeneity  can 
play  an  important  role  in  interpreting  flow,  and  especially 
transport  behavior  [see,  e.g.,  Mackay  et  al.,  1986;  Freyberg , 
1986].  Liu  et  al.  [2007]  also  showed  that  even  correct  mean 
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hydraulic  parameters  cannot  correctly  predict  drawdown 
based  on  pumping  tests  at  different  locations  in  the  domain. 
To  better  simulate  the  system  of  interest,  models  have 
grown  in  complexity  with  increasing  numbers  of  nodes  at 
a  pace  constrained  mostly  by  advances  in  computational 
power.  Often  the  number  of  nodes  in  a  model  greatly 
exceeds  the  available  number  of  measurements  resulting 
in  an  underdetermined  problem.  Methods  that  incorporate 
prior  information  or  enforce  a  structure  on  the  parameters 
result  in  an  even-determined  problem  and  the  Bayesian 
geostatistical  approach  adopted  in  this  work  is  such  a 
method. 

[3]  In  this  work,  we  illustrate  the  value  of  a  flexible 
approach  to  the  inverse  problem.  Specifically,  zonal  bound¬ 
aries  are  inferred  from  observations  assisted  by  a  small  but 
important  amount  of  prior  information.  In  the  examples  we 
evaluate,  perfect  correspondence  between  model  predictions 
and  observations  can  be  obtained  if  the  boundaries  of 
homogeneous  zones  are  known  in  advance,  but  the  real 
power  of  the  method  employed  herein  is  in  both  determin¬ 
ing  where  zonal  boundaries  are  located,  and  the  level  of 
complexity  that  should  be  present  within  each  zone.  We  also 
contend  that,  rather  than  being  a  black  box,  inverse  mod¬ 
eling  is  an  interactive  process  requiring  judgment  from  the 
practitioner  at  several  steps. 

[4]  Yeh  and  Lee  [2007]  recently  argued  that  a  paradigm 
shift  in  aquifer  characterization,  using  methods  that  allow 
for  more  detailed  parameterization,  is  essential  and  immi¬ 
nent.  We  agree  and  contend  that  for  many  studies  the 
information  available  from  wells  is  underutilized  when  used 
only  to  estimate  hydraulic  parameters  for  one  or  very  few 
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homogeneous  zones.  Li  et  al  [2007]  showed  that  hydraulic 
parameter  values  obtained  from  fitting  Theis  [1935]  curves 
to  data  from  individual  pumping  tests  at  multiple  locations 
and  interpolating  results  does  not  correctly  characterize  the 
heterogeneous  distribution  although  the  mean  parameter 
values  are  calculated  correctly.  Similarly,  Straface  et  al 
[2007]  illustrated  that  sequential  testing  with  a  distributed 
groundwater  inverse  model  often  yields  a  more  geologically 
realistic  characterization.  Inverse  modeling  allows  multiple 
measurements  of  drawdown  or  head  change  to  guide  the 
estimation  of  hydraulic  parameters.  The  inverse  problem, 
however,  is  typically  underdetermined  with  many  more 
unknown  parameters  (i.e.,  hydraulic  conductivity  values  in 
grid  cells  for  the  model)  than  measurements  (drawdown 
values  at  a  limited  number  of  wells,  for  a  limited  number  of 
pumping  tests).  Underdetermined  problems  do  not  have  a 
unique  solution.  This  shortcoming  can  be  overcome  either 
by  collecting  more  measurements  until  they  equal  or  exceed 
the  number  of  parameters  to  estimate,  by  effectively 
decreasing  the  number  of  parameters  through  lumping,  or 
by  constraining  the  solution  using  Bayes’  theorem  and  the 
inclusion  of  prior  information. 

[5]  As  Yeh  and  Lee  [2007]  suggest,  rather  than  installing 
more  wells  to  reduce  the  ratio  of  parameters  to  observations, 
multiple  pumping  tests  or  “stimulation  events”  can  be 
performed  using  various  configurations  of  existing  wells. 
The  pumping/drawdown  response  from  the  series  of  tests 
can  form  a  single  inverse  problem  that  can  be  solved  to  find 
the  set  of  hydraulic  parameters  that  best  represents  the 
observations  from  all  the  tests  simultaneously.  This  tech¬ 
nique  has  been  called  hydropulse  tomography  [To s aka  et 
al ,  1993],  a  sequential  aquifer  test  [e.g.,  Liu  et  al ,  2002],  or 
more  commonly  hydraulic  tomography  [ Bohling ,  1993; 
Gottlieb  and  Dietrich ,  1995].  Hydraulic  tomography  devel¬ 
oped  at  the  confluence  of  techniques  in  several  fields. 
Tomographic  methods  have  a  long  history  in  geophysics 
[ Sharma ,  1997]  using  seismic,  electric  current,  or  radar 
wave  propagation  through  an  aquifer  to  constrain  inversion 
of  petrophysical  properties.  By  analogy,  the  propagation  of 
water  pressure  waves  can  be  used  to  recover  hydraulic 
properties  [Tosaka  et  al ,  1993;  Vasco  etal. ,  2000;  Brauchler ; 
2005].  Several  methods  of  incorporating  this  information  into 
the  groundwater  inverse  problem  have  been  proposed, 
including  direct  inversion  of  head  data  [Tosaka  et  al , 
1993;  Gottlieb  and  Dietrich ,  1995;  Butler  et  al ,  1999;  Yeh 
and  Liu ,  2000],  steady  shape  analysis  [Bohling  et  al ,  2002], 
asymptotic  streamline  inversion  [Vasco  and  Datta-Gupta , 
1999],  and  matching  temporal  moments  of  drawdown  [Li  et 
al ,  2005;  Zhu  and  Yeh ,  2006].  Recent  work  on  both  sandbox 
models  [Illman  et  al ,  2007;  Liu  et  al ,  2007]  and  a  field 
application  [Li  et  al ,  2007]  highlight  the  power  of  hydraulic 
tomography  in  general,  and  a  Bayesian-based  parameter 
estimation  approach  specifically.  The  main  contribution  of 
the  present  work  is  in  illustrating  the  importance  that 
discontinuities  in  hydraulic  conductivity,  even  of  only  a 
single  order  of  magnitude,  can  impart  on  the  solution  of  a 
hydraulic  tomography  inverse  problem.  In  the  following 
sections,  we  discuss  methods  to  determine  and  evaluate 
zonation  candidates  to  obtain  the  best  parameter  estimation 
given  the  available  data. 

[6]  Another  motivation  for  hydraulic  tomography  is  the 
notion  that  each  pumping  test  location  in  an  aquifer  focuses 


interrogation  on  a  different  part  of  the  aquifer.  Giudici  et  al 
[1995]  examined  the  value  of  combining  multiple  stimula¬ 
tions  of  an  aquifer  to  refine  distributed  transmissivity 
estimates.  Multiple  pumping  tests  with  a  small  network  of 
wells  can  be  much  more  informative  than  a  single  pumping 
test  with  a  much  larger  network  of  observation  wells 
[Snodgrass  and  Kitanidis ,  1998].  Even  with  only  two  wells, 
Kunstmann  et  al  [1997]  showed  improved  characterization 
when  performing  unequal  strength  dipole  pumping  tests 
using  one  well  as  a  source  and  the  other  as  a  sink,  and  then 
reversing  the  configuration.  The  specific  configuration  of 
injection,  extraction,  and  monitoring  wells  in  the  field 
impacts  the  informative  coverage.  In  the  context  of  dis- 
persivity  evaluation,  Tiedeman  and  Hsieh  [2004]  explored 
the  differences  among  single-well  converging  radial  pump¬ 
ing  tests  and  equal  and  unequal  strength  dipole  tests. 
Although  their  focus  was  on  the  characterization  of  dis- 
persivity  under  different  flow  conditions,  they  also  high¬ 
lighted  the  different  areas  of  an  aquifer  interrogated  under 
the  various  flow  conditions. 

[7]  In  this  work,  we  present  a  protocol  for  implementing 
hydraulic  tomography,  and  propose  an  interactive  inversion 
scheme  using  Bayesian  geostatistical  inverse  theory  as  the 
kernel.  Prior  information  about  aquifer  parameter  character¬ 
istics  enforces  a  level  of  smoothness  or  continuity  in  the 
parameter  field.  By  “interactive”  we  mean  that  aquifer 
structure,  rather  than  being  fixed  a  priori,  is  identified 
through  the  inversion  process  and  is  evaluated  by  the 
modeler  multiple  times  both  through  the  inclusion  of  prior 
information  and  the  selection  of  zones.  This  interactive 
aspect  is  motivated  by  recognition  that  inverse  methods 
should  not  be  applied  as  a  black  box  but  rather  that 
modelers  must  be  involved  and  make  decisions  impacting 
the  results  throughout  the  process.  Intermediate  results  are 
examined  to  explore  the  impact  various  selections  of  both 
structural  parameters  (the  variogram  parameters  characteriz¬ 
ing  the  nature  of  prior  information)  and  zonation  candidates 
have  on  the  solution.  This  probabilistic  approach  imposes 
minimal  prior  assumptions  on  the  parameters,  allowing 
pertinent  information  contained  within  the  data  to  drive 
the  solution  and  explicitly  accounting  for  uncertainty. 

[8]  Uncertainty  is  introduced  both  through  epistemic 
uncertainty  and  our  inability  to  fully  characterize  inherent 
variability  of  the  media  we  model.  Epistemic  uncertainty, 
which  is  the  error  introduced  by  imperfect  and  sparse 
measurements,  roundoff  error,  approximations  of  models, 
and  conceptual  model  errors,  can  be  reduced  but  never  fully 
removed  by  improving  models  and  collecting  more  and 
better  data.  An  epistemic  error  term  accounts  for  the 
presence  of  epistemic  uncertainty  in  the  inversion.  Inherent 
variability,  or  the  heterogeneity  of  parameters  in  the  aquifer, 
cannot  be  reduced  but  is  accounted  for  through  prior 
covariance. 

[9]  The  Bayesian  geostatistical  inverse  approach  to  hy¬ 
draulic  conductivity  identification  is  usually  performed 
assuming  a  constant  but  unknown  mean  about  which  the 
best  estimate  varies  [e.g.,  Hoeksema  and  Kitanidis ,  1984; 
Snodgrass  and  Kitanidis ,  1998].  In  some  applications,  how¬ 
ever,  it  may  be  more  appropriate  to  use  a  variable  mean  if  a 
prior  trend  is  justified  [ Kitanidis ,  1997],  or  to  include  diffuse 
information  about  the  mean  [Nowak  and  Cirpka ,  2004]. 
Taking  this  concept  further,  the  known  hydrogeology  of  the 
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site  may  indicate  segregation  into  zones  with  independent 
means  that  may  improve  the  solution.  For  example,  in  the 
inversion  of  electromagnetic  borehole  flowmeter  data,  sig¬ 
nificant  improvement  was  obtained  by  assigning  zonal 
boundaries  on  the  basis  of  a  sharp  contact  between  highly 
contrasting  hydraulic  conductivity  zones  [ Fienen  et  al, 
2004].  In  some  cases,  prior  knowledge  such  as  discrete 
layering  [Bohling  et  al .,  2007],  geophysical  results 
[Hyndman  et  al ,  2000],  or  other  information  may  guide  the 
selection  of  zone  boundaries. 

[10]  It  is  critical  to  acknowledge  that  once  such  bound¬ 
aries  are  defined,  they  act  as  strong  constraints  on  the 
solution,  so  care  must  be  taken  to  justify  the  placement  of 
zones.  Ideally,  the  data  that  directly  constrain  the  inversion 
can  also  provide  enough  information  to  indicate  an  appro¬ 
priate  zonation  scheme.  Exploiting  such  information 
empowers  the  data  to  drive  the  inversion.  In  this  work, 
we  assume  only  that  a  practitioner  suspects  the  presence  of 
sharp  contacts,  and  some  expert  knowledge  is  employed  to 
select  zones  on  the  basis  of  interrogation  of  the  best  estimate 
first  obtained  assuming  a  constant  mean.  Practitioners  then 
select  among  several  candidates  for  appropriate  zonations 
and  are  guided  toward  final  selection  of  a  best  candidate 
using  metrics  based  on  examination  of  orthonormal  resid¬ 
uals  as  proposed  by  Kitanidis  [1991,  1997]. 

[n]  Given  N  distinct  screened  intervals  that  can  be  used 
either  for  pumping  or  observation,  there  are  two  main 
strategies  for  performing  the  pumping  tests  and  obtaining 
the  observations  to  be  used  for  hydraulic  tomography.  The 
first  strategy  is  to  pump  each  well  in  turn  while  monitoring 
all  the  other  wells.  The  second  strategy  is  to  perform  a 
sequence  of  dipole  tests  in  which  water  is  injected  into  one 
well  and  extracted  out  of  another  at  the  same  rate.  Again, 
the  pressure  is  observed  in  every  well  not  being  pumped. 
The  former  is  referred  to  as  “single  well  hydraulic  tomog¬ 
raphy”  (SHT)  and  the  latter  as  “dipole  hydraulic  tomogra¬ 
phy”  (DHT).  The  number  of  possible  stimulation  events  in 
the  field  for  SHT  is  N  resulting  in  TV  (N  —  1)  head 
observations.  For  DHT,  there  are  combinations  of 

dipoles  resulting  in  N(N~l)(N~2)  head  observations. 

[12]  The  number  of  stimulation  events  translates  directly 
into  the  number  of  forward  model  runs.  Thus,  for  a  given 
number  of  wells,  much  more  data  may  be  obtained  using 
DHT.  The  cost  of  installing  wells  is  greater  than  the  cost  of 
performing  pumping  tests,  so  performing  as  many  tests  as 
practical  with  a  given  well  field  should  lead  to  the  best 
characterization  possible.  For  a  given  level  of  epistemic 
uncertainty  attributed  to  each  measurement,  the  greater 
number  of  data  obtained  with  dipoles  rather  than  single 
well  tests  can  average  out  measurement  errors  thus  reducing 
uncertainty.  Furthermore,  dipoles  are  expected  to  reach 
steady  state  faster  than  single  well  tests,  so  we  explore 
DHT  in  the  synthetic  examples  of  this  work. 

[13]  The  progression  from  a  set  of  field  observations  to  a 
best  estimate  of  the  hydraulic  conductivity  field  is  guided  by 
a  combination  of  sensitivity  analysis  of  observations  to 
parameters,  the  degree  of  prior  information  weighting,  and 
the  selection  of  zones.  The  interpretation  of  the  solution  to 
the  geostatistical  inverse  problem  as  a  mapping  of  observa¬ 
tions  into  parameter  space  through  the  use  of  stochastic 
“splines”  was  discussed  by  Kitanidis  [1998]  and  further 
explored  in  the  context  of  multiple  pumping  tests  by 


Snodgrass  and  Kitanidis  [1998].  For  a  large  number  of 
wells,  it  may  be  infeasible  to  perform  all  possible  config¬ 
urations  of  pumping  tests  to  fully  realize  the  potential 
benefit  of  employing  the  DHT  protocol.  Sensitivity  matrices, 
discussed  in  section  6,  may  identify  redundant  information 
resulting  from  some  of  the  dipoles,  or  one  may  simply  select 
a  subset  of  all  dipole  configurations  to  use.  Three  synthetic 
examples  are  included:  a  stationary  heterogeneous  field  and 
two  fields  of  homogeneous  hydraulic  conductivity  that 
contain  inclusions.  The  stationary  heterogeneous  field 
serves  as  a  demonstration  model  under  ideal  conditions 
while  the  contrast  in  the  inclusion  cases  push  the  constant 
mean  assumption  to  its  limit  and  illustrates  the  value  of 
using  zonation  to  provide  for  different  mean  values  in 
distinct  zones. 

2.  Statement  of  the  Problem 

[14]  The  goal  of  this  work  is  to  perform  hydraulic 
tomography  using  several  pumping  and  observation  wells 
to  obtain  a  probabilistic  estimate  for  a  hydraulic  conductiv¬ 
ity  field.  Because  the  probabilistic  estimate  is  formulated  to 
consider  inherent  variability  and  epistemic  uncertainty,  it 
provides  a  direct  measure  of  uncertainty  about  the  most 
probable  estimate.  Using  steady  state  tests,  a  Bayesian 
geostatistical  inverse  approach  is  applied  to  the  problem. 
Prior  information  is  limited  to  the  location  of  well  screens 
available  for  pumping  and  drawdown  observations,  selec¬ 
tion  of  an  appropriate  generalized  covariance  function  used 
in  the  prior  covariance  for  the  inversion,  and  a  zonation 
assumption.  Selection  of  the  generalized  covariance  func¬ 
tion  model  is  based  on  the  general  presumed  characteristics 
(e.g.,  smoothness  or  continuity)  of  the  parameter  field.  The 
data  dictate  selection  of  the  parameters  for  the  covariance 
model  which  determines  the  degree  to  which  this  charac¬ 
teristic  is  enforced.  The  initial  zonation  assumption  is  a 
single  unknown  mean  (i.e.,  a  single  zone).  The  role  of  the 
prior  covariance  and  the  estimation  of  its  structural  param¬ 
eters  are  discussed  below.  For  this  work,  an  exponential 
generalized  covariance  function  is  used  with  a  large  integral 
scale  to  allow  the  function  to  mimic  a  linear  covariance 
function  while  behaving  as  an  appropriate  stationary 
covariance. 

[15]  The  stimulation  events  are  performed  at  steady  state 
in  two  dimensions.  The  forward  model  used  to  generate 
synthetic  observations  and  to  solve  the  inverse  problem  is 
MODFLOW-2000  (MODFLOW)  [Harbaugh  et  al ,  2000]. 
An  adjoint  state  version  of  MODFLOW  is  applied  to 
calculate  sensitivity  matrices  [ Clemo ,  2007]. 

[16]  In  steady  state  inversion,  the  aquifer  must  recover 
from  each  stimulation  event  to  ensure  that  the  observations 
from  one  test  are  not  impacted  by  conditions  still  changing 
from  the  previous  test.  Over  the  course  of  a  field  experiment 
it  may  be  difficult  to  isolate  an  aquifer  from  all  impacts  of 
regional  conditions  such  as  distant  pumping  activity,  evapo- 
transpiration,  or  recharge.  Assuming  that  regional  condi¬ 
tions  are  steady  state  over  the  long  term  can  introduce 
significant  uncertainty  and  error  in  the  estimation  of  param¬ 
eters.  Therefore  we  formulate  the  governing  equations  more 
generally  to  mitigate  changing  system  state.  If  regional 
conditions  change  at  a  timescale  slower  than  an  individual 
stimulation,  but  faster  than  the  duration  of  an  entire  hy¬ 
draulic  tomography  experiment,  this  uncertainty  can  be 
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reduced  by  working  with  drawdown  rather  than  head 
values.  By  controlling  the  conditions  under  which  data  are 
collected  the  data  value  is  increased  and  uncertainty  is 
reduced. 

[17]  Consider  the  governing  equation  of  groundwater 
flow  for  a  steady  state  isotropic  condition  [ Fetter ,  1994] 

V  •  (X(x)V0o)  =  -R(x)  (1) 

with  0O  =  <j)b  on  the  boundary  T  (2) 

where  K(x)  is  hydraulic  conductivity,  x  is  a  position  vector, 
0o  and  06  are  hydraulic  head  and  R(x)  is  the  preexistent 
recharge.  All  these  variables  are  spatially  variable. 

[is]  Next,  the  system  is  stimulated  through  pumping  in 
wells  or  any  other  stimulation,  expressed  as  a  function  of 
space  on  the  right  hand  side: 

V-0 K\7M=f(x)-R(x)  (3) 

with  0j  =  0£  on  the  boundary  T  (4) 

where  /  (x)  is  a  spatial  function  including  all  stimuli,  and  R 
(x)  remains  as  the  preexistent  recharge. 

[19]  If  we  are  interested  in  the  quantity  rf>D  =  (0 1  —  0O) 
then,  by  superposition,  we  can  subtract  the  previous  two 
equations 

V  •  (K  V(<6,))  =/(x)  -  R(x)  -  =/(x)  (5) 

with  (f)D  =  0  on  the  boundary  T  (6) 

[20]  In  this  way,  provided  that  <pb  and  R  (x)  are  constant  over 
the  timescale  of  /  (x),  we  are  left  with  a  problem  that  is 
independent  of  boundary  conditions  including  uncontrolled 
recharge  and  regional  flow.  Rather  than  measuring  and  simu¬ 
lating  head  relative  to  a  constant  datum,  the  head  measured  in 
each  well  at  the  start  of  each  stimulation  event  serves  as  the 
datum  for  that  event.  For  the  observation  data,  we  then  use 
drawdown  (head  difference)  rather  than  head.  This  is  consis¬ 
tent  with  equations  (5)  and  (6).  Previous  hydraulic  tomography 
applications  have  used  drawdown  rather  than  raw  head  data 
[see  e.g.,  Li  and  Cirpka ,  2006],  although  they  were  motivated 
by  the  need  for  integratable  data  to  calculate  moments  for 
transient  data  and  do  not  discuss  the  impact  of  regional 
conditions.  Steady  shape  conditions  described  by  Bohling  et 
al.  [2002]  provide  another  useful  alternative  intermediate 
between  transient  and  steady  state  conditions. 

[21]  To  implement  the  steady  state  method  outlined 
above,  changes  in  regional  conditions  must  occur  on  a 
timescale  slower  than  the  individual  pumping  stimulation 
event.  The  regional  conditions  may  be  different  for  each 
stimulation,  but  over  the  course  of  an  individual  stimulation, 
they  should  be  effectively  constant.  Dipoles  may  reach 
steady  state  conditions  more  rapidly. 

3.  Bayesian  Geostatistical  Inversion:  The 
Quasi-linear  Approach 

[22]  The  inverse  method  used  in  this  work  is  the  Bayesian 
geostatistical  inverse  method  [e.g.,  Hoeksema  and  Kitanidis , 


1984;  Kitanidis ,  1995].  A  summary  of  the  derivation 
follows. 

3.1.  Bayes’  Theorem 

[23]  This  method  is  rooted  in  Bayes’  theorem,  which  states 

P(*\y)  L(y\s)  p(s)  (7) 

where  p(s  |y)  is  the  posterior  probability  density  function 
(pdf)  of  the  unknown  parameters  (s),  Z(y|s)  is  the 
conditional  pdf  of  the  observations  (y)  given  the  parameters, 
referred  to  as  the  likelihood  function,  and  p( s)  is  the  prior 
pdf.  We  assume  these  pdfs  to  be  multi-Gaussian  although 
this  assumption  can  be  relaxed  using  more  computationally 
intensive  Markov  chain  Monte  Carlo  methods  [e.g., 
Michalak  and  Kitanidis ,  2003;  Fienen  et  al ,  2006].  The 
prior  pdf  contains  limited  assumptions  about  the  structure  of 
the  parameter  field  prior  to  the  experiments  being  performed 
and  measurements  being  made.  The  likelihood  function 
provides  the  means  to  update  this  prior  information  with  the 
results  of  experiments  and/or  measurements,  yielding  the 
posterior  pdf.  The  posterior  pdf  provides  the  best  estimate 
and  an  estimate  of  uncertainty  in  the  parameters  conditional 
upon  both  the  prior  information  and  the  experiments  or  data 
analysis. 

[24]  Using  a  probabilistic  approach  acknowledges  the 
nonuniqueness  of  the  parameter  estimation  problem  and 
incorporates  uncertainty  from  multiple  sources  into  the 
process  expressed  in  the  posterior  pdf.  It  is  also  possible 
to  estimate  posterior  uncertainty  of  the  estimated  parame¬ 
ters,  although  these  results  are  not  presented  in  this  work. 
Kitanidis  [1995]  describes  the  calculation  of  posterior 
covariance  necessary  to  quantify  posterior  uncertainty. 

3.1.1.  Prior  Probability  Density  Function 

[25]  The  prior  pdf  of  s  ( p(s )  in  equation  (7))  can  be 
characterized  through  its  mean  and  covariance.  We  model  s 
as  a  random  process  with  mean 

E[s]=X(3  (8) 

where  E[-]  indicates  expected  value,  0  is  a  vector  of 
unknown  drift  parameters  (p  x  1)  and  X  is  an  (m  x  p) 
matrix  of  base  functions.  The  matrix  X  associates  each  of  the 
m  values  of  s  with  the  appropriate  element  of  0,  and  can 
express  zonation  or  trend  information  known  a  priori  about  s. 
For  a  problem  with  a  single,  constant  mean,  X  is  simply  an 
(m  x  1)  vector  of  ones.  For  a  problem  with  multiple  zones, 
each  with  a  corresponding  mean  value  / 3j ,  X  is  an  assignment 
matrix  where  Xy  is  1  if  the  parameter  i  is  in  the  zone  with  the 
yth  mean  value  07-,  and  is  zero  otherwise. 

[26]  To  enforce  nonnegativity  of  s  in  unconstrained 
optimization  it  is  advantageous  to  work  with  the  logarithm 
of  hydraulic  conductivity 

s  =  In  k  (9) 

where  k  is  the  hydraulic  conductivity  vector. 

[27]  The  prior  covariance  (Q)  of  s  is 

Q(d)=E[(s-Xp)(s-X(3)T ]  (10) 

This  approach  is  an  empirical  Bayes  technique,  so  a  model 
of  Q  ( 6 )  must  be  specified  a  priori  (e.g.,  a  variogram  type), 
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but  the  parameters  (e.g.,  the  variogram  parameters,  6)  are 
estimated  on  the  basis  of  the  data.  For  the  remainder  of  this 
section,  Q(0)  is  assumed  known,  and  the  selection  of  the 
specific  model  for  Q(0)  and  its  parameters  (0)  are  discussed 
in  the  next  section. 

[28]  Assuming  a  Gaussian  distribution,  the  prior  pdf  is 


pi  s)  =  — ,  - exp 

V(27r)"det(Q) 


—  2  (s  —  X/3)rQ~‘  (s  —  X/3) 


(11) 


where  det  (•)  is  the  matrix  determinant. 

3.1.2.  Likelihood  Function 

[29]  The  likelihood  function  (X(y|s)  in  equation  (7))  is  a 
measure  of  misfit  between  predictions  of  the  measurements 
made  for  a  specific  set  of  parameter  values  weighted  by 
epistemic  uncertainty.  The  measurement  equation  relates  the 
observations  to  the  unknown  parameters 

y  =  h(s)  +  v  (12) 

where  y  is  an  (n  x  1)  vector  of  observations  (drawdown 
measurements),  s  is  an  (m  x  1 )  vector  of  unknown  parameters 
(log  hydraulic  conductivity),  h(s)  is  a  function  or  model  that 
calculates  predictions  corresponding  to  observation  values  as 
a  function  of  s,  and  v  is  an  (n  x  1)  vector  of  epistemic  error 
terms,  modeled  as  a  random  process  with  zero  mean  and 
covariance  matrix  R.  Epistemic  uncertainty  results  from 
imperfect  or  sparse  measurements  and  also  an  incomplete 
or  inappropriate  conceptual  model.  The  epistemic  error  terms 
are  assumed  independent  so  R  is  calculated  as 


[32]  A  computationally  efficient  method  to  find  the  best 
estimate  is  through  superposition  of  the  calculated  mean  and 
fluctuations 


s  =  X/3  +  QHr£ 


(17) 


where  H  is  the  sensitivity  matrix.  The  vector  of  weights  (0 


and  the  vector  of  mean  values 
cokriging  equations 


w 


are  found  through  the 


s 

HX 

y 

.  (HX)T 

0 

A 

0 

(18) 


where  Y  =  HQHr  +  R,  and  H  =  is  found  using  adjoint 
state  calculations. 

[33]  Estimating  hydraulic  conductivity  from  observations 
of  head,  h(s)  is  nonlinear  so  we  expand  the  solution  about  a 
current  best  estimate  s.  Following  the  quasi-linear  approach 
of  Kitanidis  [1995] 

h(s)  «  h(s)  +  H(s  -  s)  +  HOT  (19) 

where  HOT  are  higher-order  terms  that  are  disregarded  and 
H,  as  a  function  of  s  must  be  recalculated  at  each  lineari¬ 
zation  step. 

[34]  Updating  the  cokriging  equations  (equation  (18)), 


R  =  ^I  (13) 

where  a\  is  the  epistemic  uncertainty  parameter  and  I  is  an 
n  x  n  identity  matrix.  The  model  error  contribution  to 
epistemic  uncertainty  may,  in  reality,  be  systematic  and 
correlated  [ Gaganis  and  Smith ,  2001].  A  significant  source 
of  model  error  is  reduced  when,  as  in  this  case,  the  scale  of 
model  nodes  is  smaller  than  the  scale  of  heterogeneous 
variability  in  the  parameters  [ Gallagher  and  Doherty , 
2007].  Nonetheless,  the  decision  to  assume  uncorrelated 
homoscedacticity  in  this  work  was  made  for  practicality, 
and  the  matrix  R  can  easily  accommodate  structure  in  the 
epistemic  error  covariance. 

[30]  The  likelihood  function  is 


£(y|s) 


1 

- y  =  exp 

\J  (27t)"  det(R) 


l(y-h(s))rR  '(y-h(s)) 

(14) 


3.1.3.  Posterior  Probability  Density  Function 

[31]  Applying  Bayes’  theorem,  the  posterior  pdf  is  the 
product  of  equations  (14)  and  (11).  After  some  algebraic 
manipulation  and  removing  constants,  we  can  define  an 
objective  function  that  is  equivalent  to  the  negative  loga¬ 
rithm  of  the  posterior  pdf 

£  =  (y  —  h(s))rR-1  (y  —  h(s))  +  srGs  (15) 


E  HX 

V 

'y-h(s)  +  Hs' 

i 

o 

H 

'x 

IB, 

_ i 

A 

0 

[35]  Each  iteration  of  this  algorithm  results  in  a  new 
estimate  of  s  using  equation  (17)  with  £  and  $  calculated 
using  equation  (20)  and  an  updated  computation  of  H  For 
the  next  iteration,  we  replace  s  with  s,  and  find  a  new  s. 
Convergence  is  declared  once  ||£(s)  —  £(s)||2  or  ||s  —  s||2 
decreases  to  a  small  tolerance,  where  ||  •  H2  indicates  an  L-2 
norm. 

[36]  In  underdetermined  inverse  problems,  such  as  this 
one,  adjoint  state  calculations  dramatically  increase  the 
efficiency  of  calculating  the  sensitivity  matrix  H  which  is 
calculated  multiple  times  throughout  the  inversion  process. 
Derivations  of  the  adjoint  state  equations  are  given  by  Sykes 
et  al.  [1985]  and  Sun  [1994].  Recall  that  n  indicates  the 
number  of  observations  and  m  indicates  the  number  of 
parameters.  The  end  result  is  a  formulation  that  calculates 
H  using  n+  1  model  runs  rather  than  the  m  +  1  or  2m  +  1  runs 
required  by  the  sensitivity  equation  or  traditional  perturba¬ 
tion  sensitivities.  The  sensitivity  matrix  is  calculated  in 
physical  space  (h#  =  |jj^  but  when  we  perform  the  forward 
calculation,  we  must  convert  to  estimation  space  (^Hy  = 

in  which  the  parameters  are  log  transformed.  This  is  accom¬ 
plished  using  the  chain  rule 


where 

G  =  Q  1  -  Q_1X(XrQ_1X)_IXrQ_1  (16) 


r  _  dhi  dkj  _  dhi  9  exp  (sj) 
v  dkj  dsj  dkj  dsj 


dhi  ,  v 


(21) 


5  of  19 


W00B01 


FIENEN  ET  AL.:  BAYESIAN  HYDRAULIC  TOMOGRAPHY 


W00B01 


[37]  The  adjoint  state  calculations  are  implemented  in  a 
modified  version  of  MODFLOW.  Clemo  [2007]  provides 
the  derivation  of  the  equations  and  implementation  details 
used  for  adjoint  state  in  MODFLOW.  The  savings  in 
computational  effort  are  greater  than  tenfold  using  the 
adjoint  state  method  instead  of  finite  differences  in  the 
synthetic  cases  described  below.  Most  of  the  savings  comes 
from  the  fact  that  m  =  324  while  n  =  60,  but  additional 
overhead  savings  are  achieved  by  performing  the  adjoint 
calculations  in  a  single  call  to  the  program.  In  many 
applications,  m  is  even  greater  than  n  and  the  advantage 
to  using  adjoint  state  is  more  significant. 

3.2.  Prior  Covariance 

[38]  The  prior  covariance  in  geostatistical  inversion  is 
expressed  through  either  a  variogram  covariance  model  or  a 
generalized  covariance  function.  In  many  random  fields, 
this  covariance  model  is  adequate  to  handle  the  variability 
of  the  field,  allowing  for  clusters  of  similar  values  rising 
above  or  dipping  below  the  mean  in  a  continuously  varying 
field.  Alternatively,  sharp  discontinuities  may  indicate  zones 
of  different  mean  values  which  are  difficult  to  model  with  a 
single  covariance  and  mean  zone  assumption.  In  the  fol¬ 
lowing  subsections  we  first  discuss  the  general  covariance 
models  to  be  used  regardless  of  zonation  followed  by  a 
discussion  of  incorporation  of  zones. 

3.2.1.  Prior  Covariance  Without  Zones 

[39]  In  this  work,  we  adopt  the  exponential  covariance 
model  from  which  the  linear  covariance  model  can  be 
approximated  as  a  limiting  case 

R(h)  =  o-2  exp  (-^)  (22) 

where  a2  is  the  variance,  h  is  the  separation  distance  (in 
absolute  value)  between  nodes,  and  £  is  the  integral  scale.  A 
large  integral  scale  suggests  large  pockets  of  similar 
parameter  values,  whereas  a  small  integral  scale  indicates 
variability  at  a  small  scale  allowing  for  a  much  rougher 
solution. 

[40]  The  covariance  model  can  be  converted  to  a  variogram 

7(A)  =^(0)  ~R(h)  =  ^(l  -  exp(-*)  )  (23) 

[41]  The  generalized  covariance  function  (GCF),  is 
calculated  as 

n(h)  =  -7  (h)  +  C  (24) 

where  C  is  an  arbitrary  constant.  We  may  assume  £  is  large 
relative  to  h.  In  this  case,  variability  takes  place  over  the 
range  of  multiple  nodes.  At  the  limit  as  £  — ►  00,  exp  (—  |)  — > 
1  —  |.  Substituting  into  equation  (22) 

R(h)  =  a2  ^1  —  =  cr1  —  ^-h  =  6£  —  6h  (25) 


where  0  =  y.  Using  this  model,  since  i  >  max(/z),  the  values 
will  all  be  positive  and  the  covariance  matrix  will  be 
positive  semidefmite  and  hence  an  authorized  covariance. 
The  corresponding  variogram  is  the  linear 

7 (h)  =  R( 0)  -  R(h)  =  6£-0-0£  +  0h  =  0h  (26) 

[42]  By  assuming  £  is  constant  and  sufficiently  large  (we 
use  10  x  max  ( h )),  we  can  use  equation  (22)  with  the 
substitution  of  a2  =  6£  as  the  prior  covariance  model 

R(l ,)  =  ^exp(A)  (27) 

[43]  This  model  has  the  dual  advantages  of  having  a 
single  parameter  to  estimate  (1 9 ,  since  £  is  fixed  as  large)  and 
it  is  a  valid  covariance.  This  allows  us  to  set  lack  of 
correlation  to  zero. 

[44]  Adopting  this  GCF,  Q  in  the  prior  pdf  (p(s))  is 

Q(6»)  =  6»£exp(A)  (28) 

3.2.2.  Prior  Covariance  With  Zones 

[45]  Splitting  a  domain  into  discrete  zones  is  appropriate 
when  the  mean  values  of  various  subdomains  are  signifi¬ 
cantly  different.  Using  the  exponential  generalized  covari¬ 
ance  function  (GCF)  with  a  large  integral  scale,  sharp 
contrasts  in  conductivity  are  incompatible  with  the  smoothly 
varying  structure  described  by  the  prior  covariance.  How¬ 
ever,  sharp  contrasts  are  common  in  hydrogeologic  sedi¬ 
mentary  environments  in  the  form  of  geologic  contacts, 
particularly  at  unconformities  and  paleostream  channel 
boundaries.  As  a  result,  it  is  often  beneficial  to  segregate 
the  field  into  zones,  each  with  its  own  mean  and  uncorre¬ 
lated  from  each  other.  Determining  the  location  of  zones  is 
difficult  without  a  priori  observations  indicating  where  to 
draw  the  lines. 

[46]  The  zones  may  be  characterized  by  totally  different 
covariance  structures,  or  by  a  similar  covariance  structure 
with  different  mean  values.  Multiple  0  parameters  might 
also  be  necessary  to  characterize  anisotropy  in  which 
correlation  is  described  by  direction-specific  parameters. 
In  both  cases,  restricted  maximum  likelihood  or  another 
alternative  approach  to  finding  6  must  be  implemented.  In 
this  work,  we  assume  similar  covariance  structure  in  each 
zone  but  different  means  so  a  single  parameter  6  will  apply 
to  the  same  general  covariance  structure.  The  matrix  X 
associates  each  zone  with  its  appropriate  element  of  /?  in 
equation  (17).  The  Q  matrix,  calculated  according  to  the 
specific  GCF  used,  is  censored  such  that  for  Q y  =  0  when 
the  zth  and  yth  elements  are  in  different  zones  and  the  zones 
are  uncorrelated. 

[47]  In  some  cases,  compelling  information  regarding 
zones  is  available  in  the  form  of  boring  logs,  outcrop 
observations,  or  other  data.  Where  such  data  are  present, 
they  should  certainly  be  used  to  delineate  zones  and  reliance 
on  the  inversion  should  be  focused  on  estimating  the 
parameters  within  the  zones.  Many  approaches  are  available 
for  incorporating  outside  information  including  multiple- 
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point  geostatistics  [ Caers  and  Zhang ,  2004;  Feyen  and 
Caers ,  2006],  indicator  Kriging  [Journel,  1983],  or  cluster 
analysis  [ Tronicke  et  al .,  2004;  Bohling  et  al. ,  2007]. 
Alternatively,  there  is  often  information  in  the  data  that 
can  drive  the  selection  of  zones.  For  example,  Straface  et  al 
[2007]  used  inflections  in  observed  drawdown  curves,  as 
suggested  by  Oliver  and  Aramco  [1990]  to  indicate  bound¬ 
aries  of  different  parameter  zones.  Minimum  gradient 
support  functionals  can  also  be  used  in  a  penalization 
method  to  segregate  a  parameter  field  into  zones.  This 
approach  is  called  “image  focusing”  in  geophysics  literature 
[ Portniaguine  and  Zhdanov ,  1999,  2002]. 

3.2.3.  Thresholding  to  Select  Zones 

[48]  A  simplified  version  of  image  focusing  is  to  examine 
thresholds  within  the  best  estimate  of  a  nonzoned  model. 
Our  goal  is  to  empower  the  data  to  drive  the  inversion, 
including  the  zonation.  This  approach  allows  us  to  start  the 
inversion  clear  of  preconceived  notions  of  zonation  candi¬ 
dates  and  use  the  observations  to  indicate  appropriate  zones. 
Zones  are  indicated  by  regions  varying  about  different 
means.  On  a  histogram,  transitions  between  zones  appear 
as  gaps.  The  selection  and  definition  of  gaps  is  subject  to 
expert  interpretation,  and  several  candidates  are  selected 
interactively.  Metrics  are  derived  below  to  guide  final 
selection  among  these  candidates.  Therefore,  a  preliminary 
method  to  screen  for  zones  is  to  divide  a  histogram  of 
parameter  values  into  a  region  above  and  a  region  below  a 
threshold  identified  by  a  gap  in  the  histogram.  This  is 
illustrated  for  a  specific  example  in  Figure  la.  This  is  an 
interactive,  qualitative  approach  as  several  gaps  may  be 
present  and  expert  knowledge  is  used  to  guide  the  number 
of  zones  and  where  to  make  cutoffs.  For  example,  the  very 
definition  of  a  “gap”  on  the  histogram  is  subject  to 
interpretation.  Furthermore,  a  practitioner  may  enforce  more 
continuity  or  split  zones  by  this  method  if  soft  information 
suggests  such  actions  are  appropriate.  In  this  work,  three 
candidates  based  on  gap  locations  are  made  and  the  sol¬ 
utions  are  evaluated  for  each.  The  three  candidate  models 
were  selected  by  a  person  unfamiliar  with  the  project  and 
instructed  only  to  indicate  three  gaps  in  each  histogram.  The 
ultimate  decision  regarding  which  zonation  candidate  to 
adopt  is  based  on  cross  validation  results  (section  4)  and 
expert  discretion. 

[49]  Figure  la  shows  a  candidate  for  zonation  selected 
from  the  best  estimate  of  a  nonzoned  solution  using  the 
method  outlined  above.  Q0,  the  kernel  of  Q,  is  independent 
of  the  structural  parameter  6.  For  the  linear  generalized 
covariance  function  approximated  by  the  exponential 

Qo=fexp(A)  (29) 

[50]  To  enforce  zonation,  this  kernel  must  be  censored  by 
setting  all  covariance  matrix  values  in  Q0  to  zero  if  the  two 
nodes  are  in  different  zones.  Figure  2  shows  two  versions 
of  the  matrix  Q0:  the  top  without  zonation,  and  the  bottom 
split  into  three  zones  according  to  Figure  la.  Q  is  calcu¬ 
lated  as  0  x  Q0.  Additionally,  X  must  be  created  such  that 
nodes  are  mapped  to  the  appropriate  means  as  discussed  in 
section  3.1.1.  For  each  zonation  candidate,  the  geostatistical 


model  will  be  different,  and  the  structural  parameters,  6, 
will  need  to  be  estimated  uniquely. 


4.  Structural  Parameters 


[51]  Two  structural  parameters  must  be  estimated  along 
with  the  values  of  hydraulic  conductivity  in  the  model  nodes: 
the  epistemic  error  parameter  (oR  in  equation  (13)),  and  the 
prior  covariance  model  parameter  ( 6  in  equation  (28)). 
These  two  parameters  are  estimated  through  examination 
of  orthonormal  residuals  [ Kitanidis ,  1991,  1997].  If  the 
structural  parameters  are  appropriate,  the  drawdown  resid¬ 
uals  should  have  zero  mean  and  little  spread  about  the  mean. 
Kitanidis  [1991]  suggested  two  metrics  to  evaluate  structural 

parameters  using  orthonormal  residuals 


62 


1 


n  —  p 


n 


E 

k=P+ 1 


4 

°k 
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and 


cR  =  Q2  x  exp 


1 


n  —  p 


n 


E 

i=p+ 1 


(31) 


where  n  is  the  number  of  measurements,  p  is  the  number  of 
hydraulic  conductivity  drift  parameters,  6  are  the  residuals 
(the  actual  errors  between  observations  and  predictions  of 
drawdown),  and  a  are  the  estimation  variances.  We  seek 
parameters  that  minimize  cR  with  the  constraint  that  Q2  =  1 . 
This  is  most  easily  accomplished  by  selecting  a  value  for  o2R 
and  varying  6  over  several  orders  of  magnitude.  It  is  easier, 
a  priori,  to  determine  a  reasonable  value  for  a\  than  for  6 
since  epistemic  errors  are  related  to  measurement  errors  for 
which  there  may  be  estimates.  Otherwise,  setting  oR  to  10% 
of  the  mean  drawdown  value  is  a  reasonable  starting  value. 

[52]  The  ratio  of  controls  the  trade-off  between  fitting 

the  observations  and  enforcing  prior  information.  Small 
values  of  produce  smoother  parameter  fields.  The  optimal 

ratio  of  \  that  minimizes  cR  can  be  found  by  discrete 

aR 

graphical  intervals  as  shown  in  Figure  3  or  by  using  a 
line  search.  Q2  is  constrained  to  be  unity  by  performing 
the  minimization  of  cR  for  a  given  ratio,  then  multiplying 
both  0  and  oR  by  Q2.  At  the  next  iteration,  Q2  will  be 
equal  to  one. 

[53]  From  the  optimal  set  of  structural  parameters  for 
each  zonation  candidate,  we  can  evaluate  a  final  metric  to 
guide  selection  of  the  best  zonation  candidate.  Relative 
percent  difference  (RPD)  between  cR  and  the  optimal  oR  is 


RPD  = 


2\cR  —  cf2r 
CR  +  Or 


x  100% 


(32) 


[54]  With  the  constraint  of  Q2  =  1,  cR  is  equal  to  the 
geometric  mean  of  o2  noting 


1x“pIA,5i"'« 


n 


n 

\=p+\ 


n-p 


(33) 
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Potential  Zones  --  Bin  43 


2  4  6  8  10  12  14  16  18 

Inner  Region  Column 


Above 


Below 


Figure  1.  Zone  candidates  selected  from  an  exponential  GCF  inversion  assuming  a  single  inner  region 
zone  for  the  straight  inclusion  example.  The  left  plots  show  histograms  of  model  nodes  at  each  hydraulic 
conductivity  level.  The  right  plots  show  binary  maps  of  the  hydraulic  conductivity  fields.  The  gaps  were 
observed  at  bins  30,  43,  and  54. 


[55]  The  RPD  metric  is  therefore  a  measure  of  agreement 
between  the  epistemic  error  term  (cl)  and  the  actual  error 
calculated  from  the  orthonormal  residuals  ( cR ).  In  an 
appropriate  model,  the  epistemic  error  term  (cl)  used  in 
estimation  and  the  error  calculated  using  orthonormal  resid¬ 


uals  (c2)  should  be  close  to  each  other  since  they  are  two 
separate  measures  of  how  closely  the  measurements  should 
agree  with  the  predictions.  Selecting  the  smallest  cl  may 
lead  to  selection  of  an  overly  rough  model  that  suffers  from 
the  problem  of  overfitting.  On  the  other  hand,  RPD  is 
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Figure  2.  Comparison  of  the  Q  matrix  kernel  (Q0)  (top) 
without  zones  and  (bottom)  with  zones  for  an  exponential 
GCF  applied  to  the  straight  inclusion  synthetic  example 
with  the  lowest-threshold  zonation.  Note  that  the  lowest 
values  in  the  bottom  plot  are  equal  to  zero. 


[60]  Boundary  conditions  are  set  such  that  stimulations  in 
the  core  do  not  impact  the  boundaries.  In  this  case,  constant 
head  boundaries  were  used. 

[61 ]  Initialize  the  hydraulic  conductivity  field  as  s0.  This 
can  be  estimated  from  prior  information,  assigned  as  a 
constant  value,  or  may  be  a  previous  estimate  from  a  similar 
inversion. 

[62]  Initialize  6  and  o\.  It  is  best  to  start  with  a  high 
degree  of  prior  information  (strong  smoothing),  i.e.,  a  ratio 

of  less  than  unity,  which  makes  the  solution  smooth  and 

stable.  Smoothing  will  be  relaxed  gradually  by  increasing 
the  ratio  of  revealing  more  structure  in  the  solution  until 
the  optimal  balance  is  obtained.  This  adjustment  takes  place 
in  step  6. 

[63]  4.  Run  the  primary  and  adjoint  problem  model  runs: 
The  primary  model  run  is  a  set  of  forward  runs  (in  this  case 
using  MODFLOW)  corresponding  to  each  stimulation  given 
the  current  candidate  conductivity  field.  Initially,  the  candi¬ 
date  is  s0  and  is  subsequently  referred  to  as  s.  The  primary 
model  provides  the  current  estimate  of  the  drawdown 
values,  h(s),  which  is  a  vector  of  the  same  dimension 
as  y.  The  adjoint  problem  run  calculates  the  sensitivity 

matrix  H  =  Both  results  are  used  to  solve  equation  (15). 

[64]  5.  Find  a  new  estimate  of  s,  referred  to  as  s  by 
solving  the  following  system  of  equations: 


external  to  the  zonation  and  provides  a  measure  of  agree¬ 
ment  that  is  independent  of  the  specific  magnitude  of 
residuals.  This  provides  a  semiquantitative  guide  toward 
selection  of  an  appropriate  zonation  candidate. 


5.  Designing  a  Field  Test  and  Inversion  Protocol 

[56]  The  implementation  of  the  methodology  outlined 
above  fits  into  a  general  protocol  based  on  principles  that 
may  be  applied  to  problems  regardless  of  the  experimental 
setup  (i.e.,  SHT  or  DHT).  In  either  case,  N  wells  are 
assumed  available  which  can  each  be  configured  with  a 
pump  or  a  pressure  transducer.  This  protocol  considers 
steady  state  pumping  tests  but  transient  hydraulic  tomogra¬ 
phy  can  be  performed  with  minimal  modifications  to  the 
protocol.  The  protocol  is  outlined  as  follows: 

[57]  1.  Create  a  pumping  and  observation  protocol: 

Whether  SHT  or  DHT  will  be  performed,  a  systematic 
scheme  of  order  of  pumping  the  wells  and  collection  of 
observations  in  all  nonpumped  wells  is  required.  For  SHT, 
each  of  the  N  wells  is  pumped  individually  and  pressure  is 
observed  in  the  remaining  TV  —  1  wells.  For  DHT,  some 
subset  of  the  dipole  combinations  is  pumped  while 

pressure  is  observed  in  the  remaining  N  —  2  wells. 

[58]  2.  Perform  the  (steady  state)  pumping  tests:  Each 
pumping  test  should  be  performed  long  enough  that  steady 
state  is  achieved  at  which  point  the  drawdown  in  each 
observation  well  is  recorded.  The  datum  for  drawdown  in 
each  observation  well  is  the  head  in  that  well  immediately 
prior  to  pumping.  At  the  end  of  all  the  tests,  a  single 
combined  vector  of  drawdown  values,  y,  will  serve  as  the 
data  vector  in  equation  (12)  and  all  subsequent  calculations. 

[59]  3.  Initialize  the  forward  model: 


s  =  X/?  +  QHr£ 


(34) 


E  HX 

"y  —  h(s)  +  Hs 

- 1 

0 

H 

_ 1 

A 

0 

This  new  estimate  (s)  is  compared  with  the  previous  estimate 
(s).  These  calculations  are  iterated  until  ||£(s)  —  £(s)lk  or 
1 1 s  —  s|  1 2  decreases  beneath  a  specified  tolerance.  A  line 
search  may  be  performed  at  each  iteration  to  help  guide  the 
solution.  The  line  search  is  beneficial  for  strongly  nonlinear 
problems,  such  as  those  with  large  contrasts  [Zanini  and 
Kitanidis ,  2008].  For  many  problems,  a  line  search  is  not 


Figure  3.  Plot  of  cR  as  a  function  of  ratio  for  the  middle 
zonation  example  of  the  straight  inclusion  synthetic  case. 
The  solution  with  the  smallest  cR  value  corresponds  to  the 
optimal  trade-off  between  misfit  and  prior  information.  In 
all  cases,  Q2  is  constrained  to  unity. 
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required  unless  is  it  established  that  the  solution  is  stalled 
and  the  objective  function  does  not  decrease  in  a  satisfactory 
fashion.  An  alternative  stabilization  approach  is  the  modified 
Levenberg-Marquardt  method  [Nowak  and  Cirpka ,  2004]. 

[65]  6.  Examine  the  orthonormal  residuals  to  find  opti¬ 
mal  structural  parameters:  The  entire  process,  up  to  this 
point,  proceeds  using  the  initial  structural  parameter  values 
for  R  and  Q.  To  optimize  the  structural  parameters  the 
entire  solution  to  convergence  of  the  objective  function 
(||£(8)  —  £(s)||2)  and  the  parameters  (||s  —  s||2)  is  performed 
for  a  range  of  structural  parameters.  For  each  ratio  4-,  the 
solution  for  s  proceeds  to  convergence,  the  metrics  Q2  and  cR 
are  calculated,  and  both  structural  parameters  are  multiplied 
by  Q2,  while  maintaining  their  ratio,  so  that  the  best  estimate 
stays  the  same,  but  the  confidence  intervals  are  changed  when 
Q2  is  unity.  The  ratio  with  the  lowest  value  of  cR  is  retained  as 
the  optimal  structural  parameter  set.  When  proceeding  to  the 
next  ratio  of  4,  use  the  most  recent  s  as  the  starting  value.  It 
should  be  closer  to  the  true  s  than  a  totally  uniform  guess  and 
will  speed  convergence. 

[66]  7.  Explore  zonation  candidates:  Use  the  zonation 
approach  outlined  in  section  3.2.2. 1  above  to  select  several 
candidates  for  zonation.  For  each  candidate,  perform  the 
entire  inversion  (steps  4-6). 

[67]  8.  Evaluate  the  best  zonation  candidate:  Three  met¬ 
rics  guide  which  zonation  candidate  is  best.  The  cR  (equa¬ 
tion  (31))  value  should  be  low,  with  Q2  (equation  (30)) 
constrained  to  be  one  and  the  relative  percent  difference 
(RPD)  between  cR  and  cr|  should  be  low  (equation  (32)). 

6.  Sensitivity  Matrices  and  Splines 

[68]  Kitanidis  [1998]  and  Snodgrass  and  Kitanidis  [1998] 
discussed  the  use  of  sensitivity  matrices  and  interpolation 
splines  to  examine  the  behavior  of  geostatistical  inversion. 
The  solution  for  the  best  estimate  from  equation  (17)  at  each 
parameter  node  location  is  the  superposition  of  a  mean 
value  (X/3)  and  an  interpolation  “spline”  (QH7)  weighted 
by  £.  The  matrix  X  associates  the  node  location  with  the 
appropriate  mean  value  of  /?  which  is  estimated  by  solving 
the  system  in  equation  (18)  or  equation  (20).  The  interpo¬ 
lation  weights,  £,  are  calculated  in  the  same  system  and  are 
used,  combined  with  the  data-model  cross  covariance 
(QH  T)  to  map  the  influence  of  each  observation  in  y  to 
the  parameters.  To  understand  both  the  influence  of  specific 
measurements  in  y  on  the  best  estimate,  and  to  see  which 
nodes  are  expected  to  be  estimable,  we  can  inspect  both  the 
raw  sensitivity  (H)  and  the  splines  (QH7).  This  discussion  is 
in  the  context  of  the  specific  model  presented  in  section  7. 
The  examination  of  sensitivity  and  splines  can  guide  the 
selection  of  both  well  configurations  and  pumping  strategies. 

[69]  The  row  of  the  sensitivity  matrix  (H)  corresponding 
to  an  individual  measurement  can  be  projected  onto  the 
conductivity  field,  which  is  done  for  each  observation  for 
the  DHT  protocol  in  Figure  4.  In  Figure  4,  the  color  scale  is 
truncated  to  show  detail  of  the  field  away  from  the  well 
locations  which  would  otherwise  be  eclipsed  by  very  few 
very  high  values  immediately  adjacent  to  the  pumping  and 
observations  wells.  For  comparability,  the  reference  K  field 
is  a  homogeneous  field  with  magnitude  equal  to  the  mean  of 
the  true  solution.  A  homogeneous  reference  field  was 


selected  because  it  highlights  sensitivity  resulting  from  well 
arrangement  and  protocol  rather  than  K  contrasts. 

[70]  The  “splines”  incorporate  information  from  the  prior 
covariance  matrix  Q.  Particularly  in  the  case  of  a  linear 
prior  covariance  model,  this  information  is  diffuse  and  both 
stabilizes  the  solution  and  enforces  smoothness  in  the  best 
estimate  and  realizations.  Figure  5  shows  the  normalized 
splines  projected  onto  the  computational  domain.  The 
normalization  is  performed  relative  to  the  highest  value  in 
the  specific  column  of  QHrand  reveals  the  general  shape  of 
the  splines  throughout  the  domain.  The  splines  and  sensi¬ 
tivities  with  the  greatest  contribution  are  those  in  which  the 
injection  and  extraction  wells  span  the  domain  and  the 
observation  well  is  near  one  of  the  pumped  wells.  Examples 
are  R  4/Ob  1  and  R  12/Ob  4  in  Figure  5.  The  sign  of  both 
the  sensitivity  and  the  spline  is  controlled  by  proximity  of 
the  observation  well  to  a  pumped  well.  Negative  sensitivity 
coefficients  occur  near  extraction  wells  whereas  positive 
sensitivity  coefficients  occur  near  injection  wells.  The  least 
sensitive  configurations  are  those  when  the  dipole  is  in¬ 
duced  entirely  along  one  edge  of  the  domain  at  two  comers 
and  the  observation  well  is  either  between  or  on  the  center 
of  the  opposite  edge.  The  observations  for  Runs  2  and  14 
illustrate  these  configurations.  When  the  dipole  is  induced  at 
two  opposite  comers  with  the  observation  well  in  a  third 
comer,  sensitivity  is  also  very  diffuse  and  the  splines  not 
very  informative.  The  result  is  that  the  greatest  sensitivity 
and,  correspondingly,  the  most  informative  splines  occur 
along  the  centerline  of  the  K  field.  Interrogation  of  the 
comers  of  the  K  field  is  weaker  using  dipoles,  and  in 
general,  significant  spread  is  observed  through  the  field. 

7.  Illustrative  Examples 

[71]  We  present  three  examples  to  illustrate  the  protocol 
outlined  above.  They  are  synthetic  examples  chosen  to 
challenge  the  applicability  of  the  method,  particularly  the 
chosen  generalized  covariance  function  (GCF).  The  linear 
GCF  gives  best  results  when  applied  to  continuously 
varying  fields.  However,  it  is  a  flexible  GCF  and  often 
serves  as  the  first  choice  in  inverse  problems  when  little  is 
known  about  the  structure  of  parameter  space.  Applying  this 
model  to  a  parameter  field  with  strong  discontinuities 
highlights  its  flexibility  but  also  pushes  it  to  its  limit  of 
applicability.  Implementation  is  performed  using  a  special 
exponential  GCF  discussed  in  section  3.2.1  to  approximate 
the  linear.  Discussion  of  the  use  of  sensitivity  matrices  and 
splines  was  presented  in  the  context  of  these  specific 
examples. 

[72]  In  each  example,  the  core  is  surrounded  by  a  buffer 
zone  with  K  =  50-^-.  The  buffer  zone  extends  nearly  1  km 
in  each  direction  and  ensures  that  the  boundary  conditions 
do  not  play  a  role  in  the  tests.  In  the  stationary  heteroge¬ 
neous  example,  K  ranges  from  1  to  12^;  whereas  the 
inclusion  examples  contain  an  inclusion  with  K  =  1^ 
surrounded  by  a  K  =  10^  matrix.  The  aquifer  is  modeled 
as  confined  and  two  dimensional,  with  observations 
recorded  at  steady  state.  The  inner  region  for  ah  three 
examples  is  9  x  9  meters  with  square  discretization  and 
0.5m  nodes  resulting  in  324  K  values  to  estimate  as  shown 
in  Figure  6.  No  wells  are  included  within  the  inner  region 
but  six  wells  are  included  along  the  edges.  No  recharge  is 
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Figure  4.  Sensitivity  matrix  H  for  DHT  observations  assuming  a  homogeneous  K  field.  The  n  rows  of 
H  are  displayed  on  n  plots,  mapped  onto  the  geometry  of  the  inner  region  model  nodes.  Solid  circles  are 
extraction  wells,  solid  diamonds  are  injection  wells,  and  circled  asterisks  are  observation  wells.  The  scale 
is  truncated  at  the  high  and  low  ends  to  show  detail  in  the  middle;  actual  values  immediately  adjacent  to 
wells  are  more  extreme  than  —10  to  10.  In  each  title,  “R”  refers  to  the  dipole  test  number,  and  “Ob” 
refers  to  the  observation  number.  For  this  problem  with  six  wells,  15  dipole  tests  are  performed,  each  with 
four  observations. 
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R:  5  Ob:  1  R:  5  Ob:  2  R:  5  Ob:  3  R:  5  Ob:  4  R:  6  Ob:  1 


R:  15  Ob:  1  R:  15  Ob:  2  R:  15  Ob:  3  R:  15  Ob:  4 


Figure  5.  Interpolation  splines  for  DHT  displayed  on  the  inner  region  model  nodes.  Reference  K  field 
and  legend  are  the  same  as  in  Figure  4.  Scale  is  normalized  to  largest  value  in  each  plot  to  reveal  internal 
structure. 


applied  so  the  only  sources  of  water  is  inflow  from  the 
distant  boundaries  at  which  constant  head  is  enforced,  and 
injection  in  wells  in  the  dipole  case.  In  the  case  of  dipoles, 
the  injection  well  is  the  source  of  water  extracted  from  the 
extraction  well  and  the  head  field  is  not  perturbed  very  far 
from  the  core.  As  discussed  above,  the  changes  in  head  due 
to  pumping  (drawdowns)  are  used.  The  entire  buffer  has  a 


single  constant  value  which  is  unknown  a  priori  and  is 
estimated  along  with  the  inner  region.  As  a  result,  the  Q,  H, 
and  X  matrices  account  for  a  homogeneous  zone  that 
contains  all  nodes  in  the  buffer  zone. 

[73]  As  synthetic  examples,  the  forward  model  must  be 
run  using  the  true  K  field  to  obtain  the  measurements  (y). 
These  measurements  are  corrupted  with  a  small,  normally 
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Figure  6.  Overview  of  the  model  domain  showing  wells  and  nodes.  Each  well  can  be  either  pumped  or 
used  for  observation,  according  to  the  protocol  of  the  particular  stimulation.  The  inner  region  is 
surrounded  by  a  homogeneous  buffer  extending  1  km  in  each  direction,  depicted  in  gray,  intended  to 
attenuate  the  impact  of  boundary  conditions. 


distributed  noise  (N~(0,  1  x  10-4))  to  represent  measure¬ 
ment  uncertainty  before  being  used  in  the  inverse  solution. 
To  start  the  estimation,  the  algorithm  is  seeded  with  a 
constant  value  (the  overall  mean  of  the  K  field). 

[74]  All  examples  were  analyzed  using  both  DHT  and 
SHT  protocols  but  only  DHT  results  are  discussed  in  this 
work.  Fienen  [2007]  gives  an  expanded  discussion  of  SHT 
protocol  results.  SHT  results  were  generally  comparable 


True  Synthetic  Stochastic  Homog.  Field 


Figure  7.  True  hydraulic  conductivity  ( K)  field  for  the 
stationary  heterogeneous  field  synthetic  example. 


although  in  several  cases  the  impact  of  measurement  noise 
was  much  more  significant  than  in  the  DHT  results. 

[75]  For  the  DHT  protocol,  each  dipole  configuration  is 
used  resulting  in  15  stimulation  events,  each  with  4  obser¬ 
vations  yielding  60  discrete  drawdown  observations.  For 
each  dipole,  one  well  extracts  and  the  other  injects  at 
95  liters  per  minute.  In  dipoles  of  unequal  strength,  the 
results  of  which  well  is  injection  and  which  is  extraction 


Best  Estimate  K  Field 


Column 


Figure  8.  Best  estimate  of  the  K  field  for  the  stationary 
heterogeneous  field  synthetic  example. 
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2  4  6  8  10  12  14  16  18 


Inner  Region  Column 

Figure  9.  True  hydraulic  conductivity  ( K)  field  for  the 
straight  inclusion  synthetic  example. 

interrogate  different  areas  in  the  aquifer  and  provide  different 
information  to  the  inversion.  Therefore,  in  such  cases,  it  is 
advantageous  to  pump  each  dipole  twice;  once  for  each 
configuration.  In  this  case,  however,  the  dipoles  are  of  equal 


strength,  so  each  pair  is  pumped  only  once.  The  results  of  all 
dipoles  are  calibrated  simultaneously. 

[76]  The  protocol  outlined  above  is  then  applied  to  esti¬ 
mate  structural  parameters  and  evaluate  zonation  strategies. 

7.1.  Stationary  Heterogeneous  Field  Example 

[77]  The  first  example  is  a  smoothly  varying  synthetic 
field  shown  in  Figure  7.  This  represents  a  mildly  heteroge¬ 
neous  field  and  should  be  easy  to  find  with  this  method 
without  any  zonation.  Figure  8  shows  the  best  estimate  of 
the  hydraulic  conductivity  field.  The  result  is  smoother  than 
the  true  solution  with  the  major  features  well  represented. 
Applying  the  principle  of  parsimony,  one  could  argue  that 
an  estimation  grid  of  four  or  five  homogeneous  zones  could 
represent  this  field  nearly  as  well  as  this  highly  parameter¬ 
ized  solution.  However,  it  would  be  impossible  to  know  a 
priori  where  to  draw  the  boundaries  of  the  homogeneous 
zones.  The  flexibility  of  this  method  allowed  the  data  to 
indicate  the  appropriate  level  of  smoothness  that  could 
capture  the  variability  due  to  optimal  calculation  of  the 
structural  parameters. 

[78]  This  finding  is  possible  chiefly  because  there  are  no 
significant  discontinuities  that  the  prior  GCF  must  handle. 
This  example  proves  that  when  the  structure  conforms  to  the 
model  of  stationarity  and  spatial  continuity,  the  algorithm 
performs  very  well.  The  remaining  two  examples  are  less 
ideal  and  both  contain  significant  discontinuities  that  force 
us  to  seek  and  employ  zones  with  independent  means. 
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Best  K:  Zone  Bin  43  (middle)  Best  K:  Zone  Bin  54  (high) 


Inner  Region  Column  Inner  Region  Column 


Figure  10.  Best  estimates  for  the  straight  inclusion  synthetic  example:  (top  left)  no  zones,  (top  right) 
low-zone  candidate,  (bottom  left)  middle-zone  candidate,  and  (bottom  right)  high-zone  candidate.  The 
scales  are  different  in  each  plot  to  reveal  structure. 
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Table  1.  Structural  Parameter  Metrics  Values  for  Several  Zonation  Schemes  in  the  Straight  Inclusion  Synthetic  Example  for  DHT 


Parameter 

No  Zones 

Low  Threshold 

Middle  Threshold 

High  Threshold 

cR 

5.70  x  10"6 

1.00  x  10“6 

1.72  x  10"6 

2.33  x  10“6 

Optimal  a\ 

1.44  x  10~6 

1.19  x  10~6 

1.42  x  10~6 

1.78  x  10~6 

RPD,  % 

120 

18 

19 

27 

MSE 

19.42 

1.91  x  10“4 

9.91 

19.31 

7.2.  Straight  Inclusion  Example 

[79]  The  second  example  is  a  straight  inclusion  shown  in 
Figure  9.  Low  hydraulic  conductivity  inclusions  in  a  higher- 
conductivity  matrix  may  occur  in  a  fluvial  environment 
where  the  matrix  is  typical  river  sands  and  gravels  while  the 
lower-conductivity  inclusion  is  the  remnant  of  a  point  bar 
deposit  or  overbank  deposit  made  up  of  finer-grained 
material.  The  first  inversion  was  performed  without  using 
zonation  and  the  best  estimate  is  shown  in  the  upper  left 
panel  of  Figure  10.  The  general  pattern  is  reasonable  with 
the  boundaries  between  the  lower-^T  inclusion  and  the 
surrounding  matrix  in  the  right  location.  However,  the 
higher-^  regions  at  the  upper  and  lower  portions  of  the  inner 
region  are  artificially  smooth  and  overshoot  the  higher  K 
values  as  the  selected  GCF  is  ill  suited  to  represent  disconti¬ 
nuities.  Figure  1  shows  three  candidate  zonation  approaches 
based  on  the  first  three  significant  gaps  in  the  histogram  of  K. 
The  results  of  the  three  zonation  candidates  are  shown  in 
Figure  10.  The  lowest-zonation  candidate  in  the  upper  right 
panels  of  Figure  10  perfectly  represents  the  geometry  of  the 
true  K  regions.  The  resulting  solutions  are  nearly  perfect  as 
well  since  the  GCF  does  not  need  to  account  for  the 
discontinuity.  The  structural  parameter  metrics  are  outlined 
in  Table  1 .  The  RPD  (equation  (32))  is  lowest  for  the  low- 
zonation  candidate  in  this  case.  This  is  also  indicated  to  be 
best  by  the  mean  square  error  for  the  hydraulic  conductivity 
field  (MSE)  calculated  as 

E”!  C Kest-Ktme )2  .  . 


shown  in  Figure  13.  Remarkably,  the  angle  of  the  inclusion 
is  well  represented  in  the  solution  and  the  general  extent  is 
reasonable  with  some  smearing.  The  lowest  of  the  selected 
histogram  gaps  yields  the  closest  solution  as  indicated  by 
the  MSE  value  reported  on  Table  2.  The  associated  RPD 
(equation  (32))  values  indicate  that  zones  should  provide 
better  results  than  a  nonzoned  approach,  but  do  not  distin¬ 
guish  particularly  well  among  the  zonation  candidates. 

8.  Conclusions 

[83]  This  work  presents  an  interactive  protocol  for  hy¬ 
draulic  tomography  using  the  Bayesian  geostatistical  in¬ 
verse  method.  The  interactive  aspects  of  the  protocol  allow 
the  imposition  of  expert  judgment  and  allow  practitioners  to 
investigate  the  ramifications  of  their  choices.  The  metrics 
from  orthonormal  residuals  provide  important  guidance 
regarding  selection  of  zones  and  smoothing  versus  data 
reproduction.  However,  it  is  up  to  practitioners  to  evaluate 
the  geologic  realism  provided  in  each  solution.  In  the 
determination  of  structural  parameters  and,  to  a  greater 
extent,  in  the  selection  of  zones,  interactive  opportunities 
allow  practitioners  to  audit  to  algorithm  and  override 
“decisions”  indicated  by  the  metrics  available.  Inverse 
modeling  is  not  a  black  box,  and  expert  knowledge  still 
plays  an  important  role  in  obtaining  realistic  and  meaningful 
estimates  of  hydrogeologic  parameters. 

[84]  The  Bayesian  method  explicitly  considers  epistemic 
uncertainty  and  imposes  a  minimum  of  prior  assumptions 
on  the  parameters.  We  showed  that  favorable  results  can  be 


where  Kest  and  Ktrue  are  the  estimated  and  actual  parameter 
values  in  each  of  the  m  model  nodes  of  the  inner  region. 

[so]  This  calculation  is  only  possible  for  synthetic  cases, 
of  course.  The  ability  to  compare  the  results  among  zona¬ 
tion  candidates  to  “truth”  is  an  advantage  of  evaluating  the 
method  using  a  synthetic  example. 

7.3.  Angled  Inclusion  Example 

[si]  The  third  example,  depicted  in  Figure  1 1 ,  is  a  narrow 
inclusion  at  an  angle  oblique  to  the  principal  directions  of 
the  model  domain  which  could  indicate  a  clay  lens.  Three 
zonation  candidates  were  chosen  on  the  basis  of  the  histo¬ 
grams  presented  in  Figure  12.  The  gaps  are  not  as  clearly 
defined  as  in  the  straight  inclusion  example.  The  lowest  gap 
was  initially  selected,  but  the  inclusion  was  implied  to  be 
too  small  and  the  imposition  of  such  an  erroneous  zonation 
scheme  was  unstable  precluding  a  meaningful  inversion. 
The  process  of  selection  and  evaluation  of  the  zones  is 
subject  not  only  to  the  metrics  presented  in  this  work,  but 
also  judgment  of  the  practitioner. 

[82]  As  with  the  previous  example,  the  general  shape  of 
the  inclusion  is  found  with  DHT  even  without  zones  as 
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Figure  11.  True  hydraulic  conductivity  ( K)  field  for  the 
angled  inclusion  synthetic  example. 
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Figure  12.  Zone  candidates  selected  from  an  exponential  GCF  inversion  assuming  a  single  inner  region 
zone  for  the  angled  inclusion  example.  The  left  plots  show  histograms  of  pixels  at  each  hydraulic 
conductivity  level.  The  right  plots  show  binary  maps  of  the  hydraulic  conductivity  fields.  The  gaps  were 
observed  at  bins  18,  23,  and  34. 


obtained  with  a  relatively  small  number  of  wells,  even  if 
located  entirely  outside  the  domain  in  which  hydraulic 
conductivity  parameters  are  estimated.  This  could  be  the 
field  setup  in  a  case  of  limited  access  when  seeking  lateral 
extent  information,  or  a  series  of  cross-sectional  multilevel 


screened  wells  in  a  cross-hole  context  analogous  to  seismic 
or  electrical  resistivity  tomography.  In  the  cross-sectional 
case,  a  three-dimensional  model  may  be  more  appropriate. 
However,  our  method  is  not  limited  to  cases  where  all  wells 
are  outside  the  domain  of  interest.  In  many  cases,  wells  are 
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Best  K:  Zone  Bin  23  (middle)  Best  K:  Zone  Bin  34  (high) 
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Figure  13.  Best  estimates  for  the  angled  inclusion  synthetic  example:  (top  left)  no  zones,  (top  right) 
low-zone  candidate,  (bottom  left)  middle-zone  candidate,  and  (bottom  right)  high-zone  candidate.  The 
scales  are  different  in  each  plot  to  reveal  structure  within  each  solution. 


available  within  the  domain  of  interest  and  these  wells  may 
be  extremely  valuable,  providing  significant  information 
about  the  parameters.  An  important  source  of  uncertainty 
is  alleviated  in  this  protocol  by  referencing  all  head  meas¬ 
urements  to  immediately  antecedent  levels.  Using  draw¬ 
down  mitigates  the  impact  that  longer-term  variations  in 
regional  conditions  could  otherwise  have  on  confounding 
the  parameter  estimation  over  the  relatively  short  duration 
of  each  hydraulic  tomography  stimulation  event. 

[85]  The  adjoint  state  implementation  is  an  efficient  and 
accurate  means  to  calculate  the  sensitivities.  Examination  of 
sensitivities  and  interpolation  splines  provides  valuable 
information  about  which  well  configurations  provide  the 
most  information  about  the  parameter  field.  Particularly  in 
the  case  of  dipole  hydraulic  tomography,  this  information 
can  guide  the  selection  of  a  subset  of  dipoles  from  an 
existing  well  field  to  incorporate  into  a  hydraulic  tomogra¬ 
phy  experiment. 

[86]  The  examination  of  orthonormal  residuals  proved  to 
be  an  effective  tool  in  evaluating  the  trade-off  between  data 


reproduction  and  prior  information.  The  cR  and  Q2  metrics 
were  supplemented  by  a  new  metric,  the  RPD,  which  is  the 
relative  percent  difference  between  the  optimal  epistemic 
error  parameter  (cr|)  and  cR,  interpreted  as  a  proxy  for  a\ 
measured  from  a  given  solution.  The  RPD  should  be  low 
when  these  values  are  in  agreement  and  it  was  shown  to 
provide  guidance  about  when  to  accept  candidate  models 
(on  the  basis  of  zonation  choices)  to  improve  the  best 
estimate  results. 

[87]  Zonation  was  shown  to  be  a  powerful  force  in 
obtaining  best  estimates  for  problems  of  low-if  inclusions. 
The  best  estimate  of  the  K  field  obtained  without  zones 
resulted  in  oversmoothing  where  the  prior  covariance  model 
was  taxed  too  much  trying  to  capture  discontinuities  in  the 
field.  Zones  determined  through  active  examination  of 
preliminary  nonzoned  results  yielded  excellent  solutions, 
particularly  when  the  zone  choices  closely  matched  the 
inclusion  dimensions.  When  the  zones  were  selected  inap¬ 
propriately,  the  solution  suffered,  but  the  metrics  used  here 
provided  guidance  on  which  to  accept  and  which  to  reject. 


Table  2.  Structural  Parameter  Metrics  Values  for  Several  Zonation  Schemes  in  the  Angled  Inclusion  Synthetic  Example  for  DHT 


Parameter 

No  Zones 

Low  Threshold 

Middle  Threshold 

High  Threshold 

cR 

7.66  x  1(T6 

4.38  x  10~6 

4.37  x  10~6 

4.16  x  10~6 

Optimal  (Jr 

1.24  x  10"6 

1.38  x  10“6 

1.50  x  10"6 

1.32  x  10“6 

RPD,  % 

144 

104 

98 

104 

MSE 

60.77 

9.44 

10.84 

19.59 
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The  promise  shown  by  a  combination  of  the  linear  prior 
covariance  model  and  zones  motivates  future  development 
of  more  robust  zonation  strategies. 
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